A micromechanics-based analytical solution for the effective thermal conductivity of composites with orthotropic matrices and interfacial thermal resistance

Abstract

We obtained an analytical solution for the effective thermal conductivity of composites composed of orthotropic matrices and spherical inhomogeneities with interfacial thermal resistance using a micromechanics-based homogenization. We derived the closed form of a modified Eshelby tensor as a function of the interfacial thermal resistance. We then predicted the heat flux of a single inhomogeneity in the infinite media based on the modified Eshelby tensor, which was validated against the numerical results obtained from the finite element analysis. Based on the modified Eshelby tensor and the localization tensor accounting for the interfacial resistance, we derived an analytical expression for the effective thermal conductivity tensor for the composites by a mean-field approach called the Mori-Tanaka method. Our analytical prediction matched very well with the effective thermal conductivity obtained from finite element analysis with up to 10% inhomogeneity volume fraction.

Introduction

Due to their exceptional effective thermal conductivity, reinforced composites are widely used in various applications, such as electronic devices^1,2, thermal energy storage^3,4, and thermoelectric devices⁵. For example, the thermal conductivities of epoxy reinforced with silicon⁶, graphite^7,8, or aluminium nitride⁹ are significantly improved compared to pure epoxy. These highly conductive composites have been investigated for use as efficient thermal heat sinks in the electronic packaging design^10,11,12. Meanwhile, composites made of low conductive materials can serve as good insulators^13,14. To efficiently design and use these composites, it is essential to understand and be able to predict the effective thermal conductivity of the composites.

Numerous previous studies have focused on the effective material properties of composites based on numerical simulations and theories^{15,16,17,18,19,20,21,22,23,24,25,26}. The effective thermal conductivities of macroscale composites have been computed by finite element analysis (FEA) as a function of the shape and orientation of their inhomogeneities^15,16,27. Molecular dynamics simulations were used to predict the effective thermal conductivities of composites possessing nanoscale inhomogeneities^23,26. However, multiple calculations are required to obtain statistically meaningful values of effective conductivity by changing the positions, numbers, and the orientation distributions of the inhomogeneities in a simulation cell^16,28. For example, when the thermal conductivity of a spherical inhomogeneity is twenty time larger than that of a matrix, an FEA study based on a representative volume element (RVE) containing ten inhomogeneities showed that the standard deviation of thermal conductivity became as large as 10% of the predicted thermal conductivity¹⁶. Sufficiently large RVE calculations result in smaller statistical errors, but they consume significant computational resources.

Alternatively, homogenization theories have been applied to predict the effective thermal conductivity of reinforced composites with a relatively low (<20%) volume fraction of inhomogeneities^{4,21,22,24,25}. In a micromechanics-based approaches to elasticity, the effective stiffness of a composite can be computed when the Eshelby tensor of a given inclusion is known²⁹. The Eshelby tensor relates eigenstrain and constrained strain in the single inclusion problem, and can be applied to solve inhomogeneity problems by considering the equivalent eigenstrain. The Eshelby tensor has been theoretically calculated for a variety of inclusion shapes, including prolate spheroids, oblate spheroids, spheres, and cylinders in the case of an isotropic matrix³⁰. The eigen-intensity problem in heat conduction is mathematically analogous to the eigenstrain problem in elasticity^22,30 (See Supplementary Table 1). Hence, the Eshelby tensor is constant for the ellipsoidal inclusion in heat conduction in the absence of an interfacial resistance, as it is constant for analogous elasticity problems. The effective thermal conductivity of the composites with multiple inhomogeneities can be expressed in terms of the Eshelby tensor and the localization tensor (analogous to the stress concentration tensor in elasticity), based on homogenization approaches such as Mori-Tanaka or self-consistent methods^21,22,24,25.

Theoretical studies for the isotropic matrix and ellipsoidal inclusion without interfacial thermal resistance have well been established^19,21,30. However, in realistic composites, interfacial thermal resistance arises from many reasons, such as relative roughness, lattice mismatch, and poor chemical or mechanical adhesion^31,32, which results in a temperature jump across the interface. Such interfacial thermal resistance, the so-called Kapitza resistance, reduces the thermal conductivity of the composites^17,20,33. Another complexity arises because the thermal conductivity of a matrix is anisotropic, i.e., thermal conductivity changes with the direction of the heat flux applied to matrix. For example, single-crystal metals and ceramics in a tetragonal or hexagonal crystal lattice have transversally isotropic thermal conductivity where the conductivities along two directions are identical but different to the conductivity along the other direction^34,35,36. Polymer matrices synthesized through extrusion and drawing processes can also have transversally isotropic conductivity because the polymer chains are aligned along one direction³⁷. They have different axial and lateral physical properties depending on processing conditions, such as the processing temperature and extrusion rate^37,38,39. Orthorhombic crystals such as cementite⁴⁰, titanium alloy⁴¹, and tin selenide (SnSe)⁴² have three independent material constants in the thermal conductivity tensor⁴³. Crystal families with less symmetry can have non-zero off-diagonal components in the thermal conductivity tensor, i.e., the direction of heat flow may not be exactly same as the direction of temperature gradient. There are several experimental studies on the thermal conductivity of composite materials with anisotropic matrices^44,45. However, existing theoretical studies consider either anisotropic matrices with zero interfacial resistance⁴⁶, or isotropic matrices with finite interfacial resistance²². To the best of the authors’ knowledge, there exists no theoretical study simultaneously considering the anisotropy of matrix and the interfacial thermal resistance.

In this work, we derive, for the first time, an analytical expression for the effective thermal conductivity of composites composed of spherical inhomogeneities and orthotropic matrices by accounting for the effect of interfacial thermal resistance. In the first half of the remainder of this paper, we obtain the analytical expression of the Eshelby tensor for an eigen-intensity problem when the matrix is an orthotropic material in the absence of the interfacial resistance. Based on the Eshelby tensor, we compute the heat flux within a single inhomogeneity and the effective modulus of composite. In the second half, we obtain a modified Eshelby tensor that accounts for the interfacial thermal resistance. For the single inhomogeneity problem, theoretical predictions of the heat flux within the inhomogeneity under external heat flux, and the amount of the temperature jump at the interface, are validated against numerical calculations based on FEA for a wide range of interfacial thermal resistances. We then apply a micromechanics-based homogenization method to derive an analytical solution of the effective thermal conductivities along three axes for the particle-reinforced composites with orthotropic matrices. The effective thermal conductivity prediction correctly converges to that of a porous matrix at the infinite interfacial resistance limit and that of perfect interface solution in the zero interfacial resistance limit. We show that our analytical prediction matches very well with the FEA results for an RVE of particle-reinforced composites.

Results and Discussion

Effective Thermal Conductivity in the Absence of Interfacial Resistance

In the steady state, the governing equation for heat conduction is written as,

$${K}_{{0}_{ij}}\frac{\partial }{\partial {x}_{i}}(\frac{\partial T({\boldsymbol{x}})}{\partial {x}_{j}})+g({\boldsymbol{x}})=0$$

(1)

where g is a heat source, T is temperature field, and K₀ is a symmetric second order thermal conductivity tensor. The repeated small letter indices represent the dummy indices that imply summation over all the values from 1 to 3. The Green’s function G(x − y) in the steady state heat conduction equation is defined as the temperature field at a position x in the presence of a unit heat source at another position y in an infinite medium,

$${K}_{{0}_{ij}}\frac{{\partial }^{2}G({\boldsymbol{x}}-{\boldsymbol{y}})}{\partial {x}_{i}\partial {x}_{j}}+\delta ({\boldsymbol{x}}-{\boldsymbol{y}})=0.$$

(2)

By solving the equation, the Green’s function is obtained as,

$$G({\boldsymbol{x}}-{\boldsymbol{y}})=\frac{1}{4\pi \sqrt{{\rm{\det }}({{\boldsymbol{K}}}_{0})[{({\boldsymbol{x}}-{\boldsymbol{y}})}^{T}{{\boldsymbol{K}}}_{0}^{-1}({\boldsymbol{x}}-{\boldsymbol{y}})]}}.$$

(3)

(See Supplementary Note 1 for the details). For the isotropic materials $({K}_{ij}={k}_{0}{\delta }_{ij})$, the Greens’ function reduces to a well-known function that is inversely proportional to the distance, $\frac{1}{4\pi {k}_{0}|{\boldsymbol{x}}-{\boldsymbol{y}}|}.$ In the absence of interfacial thermal resistance, the Green’s function is introduced to derive the Eshelby tensor S_ik of an eigen-intensity problem that relates the intensity field ${\boldsymbol{e}}=-\,\nabla T$ and the eigen-intensity field e^* within the inclusion as^30,47,

$${S}_{ik}({\boldsymbol{x}})=\frac{\partial }{\partial {x}_{i}}{\int }_{V}\frac{\partial G({\boldsymbol{x}}-{\boldsymbol{y}})}{\partial {y}_{j}}d{\boldsymbol{y}}\,{K}_{{0}_{jk}}\,{\rm{and}}\,{e}_{j}={S}_{ij}{e}_{j}^{\ast }$$

(4)

where V is the volume of an inclusion whose thermal conductivity is identical to the matrix (see Fig. 1(a)). The eigen-intensity field e^* can be considered as a fictitious temperature gradient field produced without external heat flux for an isolated inclusion, and the intensity field e is the temperature gradient within the inclusion when it is embedded in the matrix. With the Eshelby tensor, we can solve the eigen heat flux problem, which relates the heat flux within the inclusion q and eigen heat flux q^* as,

$${q}_{i}={C}_{ij}{q}_{j}^{\ast }.$$

(5)

where the C is known as the conjugate Eshelby tensor,

$${\boldsymbol{C}}=-\,{{\boldsymbol{K}}}_{0}{\boldsymbol{S}}{{\boldsymbol{K}}}_{0}^{-1}+{\boldsymbol{I}}.$$

(6)

By adopting the classical potential theory and the mathematical analogy with electrostatics⁴⁸ (See Appendix B in the reference), one can simplify Eq. (4) for an ellipsoidal inclusion as

$${\boldsymbol{S}}=\frac{{\rm{\det }}({\boldsymbol{a}})}{2}{\int }_{0}^{\infty }\frac{{({{\boldsymbol{a}}}^{2}+s{{\boldsymbol{K}}}_{0})}^{-1}{{\boldsymbol{K}}}_{0}}{\sqrt{\det ({{\boldsymbol{a}}}^{2}+s{{\boldsymbol{K}}}_{0})}}ds.$$

(7)

The tensor ${\boldsymbol{a}}=[\begin{array}{ccc}{a}_{1} & 0 & 0\\ 0 & {a}_{2} & 0\\ 0 & 0 & {a}_{3}\end{array}]$ is a diagonal matrix with its diagonal elements being the half the length of the principal axes of an ellipsoidal inclusion whose volume is defined as $\{({x}_{1},{x}_{2},{x}_{3}):{\sum }_{i=1}^{3}\frac{{x}_{i}^{2}}{{a}_{i}^{2}}\le 1\}$. For example, for the spherical inclusion, a can be expressed as ${a}_{ij}=R{\delta }_{ij}$, where R is the radius of inclusion.

Figure 1

Schematic of (a) an eigen-intensity problem, (b) a single inhomogeneity problem and (c) interfacial thermal resistance. (d) The temperature jump at the interface.

In this study, we consider the Eshelby tensor for a spherical inclusion in the orthotropic matrix with three independent thermal conductivity coefficients. The thermal conductivity tensors of matrix (K₀) and inhomogeneity (K₁) are defined as follows,

$${{\boldsymbol{K}}}_{0}=[\begin{array}{ccc}{k}_{1} & 0 & 0\\ 0 & {k}_{2} & 0\\ 0 & 0 & {k}_{3}\end{array}]\,{\rm{a}}{\rm{n}}{\rm{d}}\,{{\boldsymbol{K}}}_{1}=[\begin{array}{ccc}{\kappa }_{1} & {\kappa }_{12} & {\kappa }_{13}\\ {\kappa }_{12} & {\kappa }_{2} & {\kappa }_{23}\\ {\kappa }_{13} & {\kappa }_{23} & {\kappa }_{3}\end{array}].$$

(8)

Now, we derive a closed form expression of S_ij for a spherical inclusion by plugging the thermal conductivity tensor K₀ into the Eq. (7) with ${a}_{ij}=R{\delta }_{ij}$ as,

$${S}_{IJ}={\delta }_{IJ}\frac{{R}^{3}}{2}{\int }_{0}^{\infty }\frac{{({R}^{2}+s{k}_{I})}^{-1}{k}_{I}}{\sqrt{{\prod }_{\ell =1}^{3}({R}^{2}+s{k}_{\ell })}}ds={\delta }_{IJ}\frac{{({k}_{I})}^{3/2}}{2\sqrt{{\rm{\det }}({{\boldsymbol{K}}}_{0})}}{\int }_{0}^{\infty }\frac{1}{(s^{\prime} +1)\sqrt{{\prod }_{\ell =1}^{3}(s\text{'}+\frac{{k}_{I}}{{k}_{\ell }})}}ds\text{'}$$

(9)

where for the last equality, we set $s^{\prime} =s{k}_{I}/{R}^{2}$. Because both k_I and ${R}^{2} > 0$, range of integration does not change (See Supplementary Note 2 for the details). The repeated capital indices I and J, are not summed over. We note that S is a diagonal matrix because K₀ is a diagonal matrix. By defining the anisotropy factors of the matrix as $A={k}_{I}/{k}_{L},\,B={k}_{I}/{k}_{M}$ with $I\ne L,I\ne M,$ and $L\ne M$, we can simplify the Eshelby tensor as

$${S}_{IJ}={\delta }_{IJ}\frac{{({k}_{I})}^{\frac{3}{2}}}{\sqrt{{\rm{\det }}({{\boldsymbol{K}}}_{0})}}\frac{1}{A-B}(\frac{E({\cos }^{-1}(\sqrt{B}),\frac{A-1}{B-1})}{\sqrt{1-B}}-\frac{E({\cos }^{-1}(\sqrt{A}),\frac{B-1}{A-1})}{\sqrt{1-A}})$$

(10)

where $E(\theta ,m)={\int }_{0}^{\theta }\sqrt{1-m\,{\sin }^{2}\theta ^{\prime} }d\theta ^{\prime}$ is elliptic integral of the 2^nd kind (See Supplementary Note 3 for the details). For example, when $I=J=2$, either $L=1,\,M=3$ or $L=3,\,M=1$ is used for calculating the S₂₂ component, and either L, M combination results in the same result. We note that, at the limit $B\to A\,({\rm{or}}\,A\to B)$, the Eshelby tensor reduces to the transversely isotropic matrix result in a closed form solution (See Supplementary Note 3). At the limit $A,B\to 1$, the Eshelby tensor converges to an isotropic matrix case such that ${S}_{ij}=\frac{1}{3}{\delta }_{ij}$ ^22,30_, (see Supplementary Fig. S1). The three independent values (S₁₁, S₂₂, S₃₃) are plotted in terms of ${k}_{1}/{k}_{2}$ and ${k}_{1}/{k}_{3}$, and we validate our analytical solutions against the numerical evaluation of Eq. (4) (see Supplementary Fig. S1).

It has been proven that the heat flux ${{\boldsymbol{q}}}^{{\rm{inh}}}$ within an ellipsoidal inhomogeneity embedded in an infinite matrix under the presence of a constant external far field heat flux ${{\boldsymbol{q}}}^{{\rm{ext}}}$ is uniform, and so is the Eshelby tensor S, regardless of the materials symmetry of the matrix^29,30 (see Fig. 1(b)). Heat fluxes, ${{\boldsymbol{q}}}^{{\rm{inh}}}$ and ${{\boldsymbol{q}}}^{{\rm{ext}}}$, are related by the localization tensor B^22,30 as $\,{{\boldsymbol{q}}}^{{\rm{inh}}}={\boldsymbol{B}}{{\boldsymbol{q}}}^{{\rm{ext}}}$, where ${\boldsymbol{B}}={[{\boldsymbol{I}}+{\boldsymbol{C}}{{\boldsymbol{K}}}_{0}({{\boldsymbol{K}}}_{1}^{-1}-{{\boldsymbol{K}}}_{0}^{-1})]}^{-1}$. Using the Eshelby tensor in the orthotropic matrix in Eq. (10), we predict the heat flux within the inhomogeneity with any arbitrary thermal conductivity tensor K₁. For an inhomogeneity with isotropic or cubic symmetry whose thermal conductivity tensor is given as ${K}_{{1}_{ij}}=\kappa {\delta }_{ij}$, the heat flux expression can be simplified as

$${q}_{I}^{{\rm{i}}{\rm{n}}{\rm{h}}}=\frac{\kappa }{{k}_{I}+{S}_{II}(\kappa -{k}_{I})}{q}_{I}^{{\rm{e}}{\rm{x}}{\rm{t}}}$$

(11)

Our solution is validated against the numerical calculations based on FEA where we consider a single inhomogeneity surrounded by an orthotropic medium in a cubic shape. We set the edge of the medium to be 15 times larger than the diameter of the inhomogeneity to reasonably describe the infinite medium²² as depicted in Fig. 2. The FEA is performed by COMSOL software with a total of 416,606 tetrahedral quadratic elements in the matrix and inhomogeneity. In this simulation, the unit heat flux boundary conditions are considered in the x, y, and z directions to study the effect of the anisotropy (${q}_{I}^{{\rm{ext}}}=1\,{\rm{W}}/{{\rm{m}}}^{2},{q}_{J}^{{\rm{ext}}}=0$ with $I\ne J$). We carry out the calculations using ${k}_{1}=1,{k}_{2}=2,{k}_{3}=3$ and $\kappa =10\,({\rm{W}}/{\rm{mK}})$. As expected, the calculated heat flux within the inhomogeneity is uniform and dependent on the external heat flux direction (see Fig. 2(b)), and matches very well with the theoretical prediction (see Fig. 3).

Figure 2

(a) The mesh configuration for a single inhomogeneity and the matrix used in FEA. (b) Heat flux distribution within the matrix and the inhomogeneity at x, y, zplane for three different heat flux directions. The thermal conductivities used for the results are ${k}_{1}=1,{k}_{2}=2,{k}_{3}=3$and $\kappa =10$ (W/mK).

Figure 3

Normalized uniform heat flux values ((a) q₁/q₀, (b) q₂/q₀, (c) q₃/q₀) within the inhomogeneity as a function of anisotropy factors where q₀ is the magnitude of the heat flux at the boundary. The isotropic thermal conductivity of a single inhomogeneity is 10 (W/mK).

The effective thermal conductivity of a composite with multiple inhomogeneities can be predicted by considering the interaction between the inhomogeneities. In a mean field approach, as in the Mori-Tanaka method, the heat flux within the inhomogeneity is related to the average heat flux within the matrix. The Mori-Tanaka model is known to predict effective properties well at a relatively low inhomogeneity volume concentration (<20%) and is more convenient than the self-consistent method, which relies on a nonlinear implicit equation. In the absence of the interfacial resistance, the effective thermal conductivity based on the Mori-Tanaka method can be obtained as^19,30

$${K}_{I}^{{\rm{eff}}}=\frac{{k}_{I}[(1-{S}_{II}){c}_{0}{k}_{I}+({c}_{1}+{S}_{II}{c}_{0})\kappa ]}{[{k}_{I}+{S}_{II}{c}_{0}(\kappa -{k}_{I})]}$$

(12)

here, c₀ and c₁ refer to volume concentrations of the matrix and inhomogeneity, respectively; thus, ${c}_{0}+{c}_{1}=1$.

Effective Thermal Conductivity in the Presence of Interfacial Resistance

We now turn our attention to the realistic system, where interfacial thermal resistance is present^31,32. The interfacial thermal resistance α is defined as

$${T}^{{\rm{out}}}-{T}^{{\rm{in}}}=-\,\alpha {\boldsymbol{q}}\cdot {\boldsymbol{n}}$$

(13)

where ${T}^{{\rm{out}}}$ and ${T}^{{\rm{in}}}$ refer to temperatures at the outer and inner surfaces of the interface, respectively, q is the heat flux at the interface, and the n is the outward surface normal vector (see Fig. 1(c,d)). The SI unit of interfacial thermal resistance α is [m²K/W].

The interfacial resistance augments an additional surface integration term in the eigen-intensity problem, as follows,

$${e}_{m}({\boldsymbol{x}})=\frac{\partial }{\partial {x}_{m}}{\int }_{V}\frac{\partial G({\boldsymbol{x}}-{\boldsymbol{y}})}{\partial {y}_{i}}d{\boldsymbol{y}}{K}_{{0}_{ij}}{e}_{j}^{\ast }+\frac{\partial }{\partial {x}_{m}}\alpha {K}_{{0}_{ij}}{K}_{{0}_{sr}}{\int }_{\partial V}\frac{\partial G({\boldsymbol{x}}-{\boldsymbol{y}})}{\partial {y}_{j}}{n}_{i}({\boldsymbol{y}}){n}_{s}({\boldsymbol{y}})({e}_{r}({\boldsymbol{y}})-{e}_{r}^{\ast }({\boldsymbol{y}}))d{\boldsymbol{y}}.$$

(14)

It has been found that the heat flux within a spherical inclusion is uniform in the presence of interfacial thermal resistance²². Although a previous study²² claims that the heat flux within an elliptical inclusion is also uniform in the presence of interfacial thermal resistance, our numerical tests reveal that it is non-uniform (See Supplementary Note 4). Because the intensity field in the spherical inclusion is uniform, we can simplify Eq. (14) as follows:

$$\begin{array}{rcl}{e}_{m}({\boldsymbol{x}}) & = & \frac{\partial }{\partial {x}_{m}}{\int }_{V}\frac{\partial G({\boldsymbol{y}}-{\boldsymbol{x}})}{\partial {y}_{i}}{\rm{d}}{\boldsymbol{y}}{K}_{{0}_{ij}}{e}_{j}^{\ast }\\ & & +\alpha {K}_{{0}_{ij}}{K}_{{0}_{sr}}\frac{\partial }{\partial {x}_{m}}{\int }_{\partial V}\frac{\partial G({\boldsymbol{y}}-{\boldsymbol{x}})}{\partial {y}_{j}}{n}_{i}({\boldsymbol{y}}){n}_{s}({\boldsymbol{y}}){\rm{d}}{\boldsymbol{y}}({e}_{r}-{e}_{r}^{\ast })\\ & = & {S}_{mj}{e}_{j}^{\ast }+\alpha {K}_{{0}_{ij}}{K}_{{0}_{sr}}{M}_{ijms}({e}_{r}-{e}_{r}^{\ast })\end{array}$$

(15)

$\text{where}\,$

$${M}_{ijms}\equiv \frac{{\rm{\partial }}}{{\rm{\partial }}{x}_{m}}{\int }_{{\rm{\partial }}V}\frac{{\rm{\partial }}G({\boldsymbol{y}}-{\boldsymbol{x}})}{{\rm{\partial }}{y}_{j}}{n}_{i}({\boldsymbol{y}}){n}_{s}({\boldsymbol{y}}){\rm{d}}{\boldsymbol{y}}.$$

From ${e}_{m}={S}_{mj}^{M}{e}_{j}^{\ast }$, the modified Eshelby tensor is given as

$${{\boldsymbol{S}}}^{M}={({\boldsymbol{I}}-\alpha {{\boldsymbol{K}}}_{0}:{\boldsymbol{M}}{{\boldsymbol{K}}}_{0})}^{-1}({\boldsymbol{S}}-\alpha {{\boldsymbol{K}}}_{0}:{\boldsymbol{M}}{{\boldsymbol{K}}}_{0}).$$

(16)

Where the colon in Eq. (16) is a double dot product. We also find that the Eq. (16) can be further simplified (See Supplementary Note 5 for details), as follows:

$${{\boldsymbol{S}}}^{M}={({\boldsymbol{I}}+\frac{\alpha }{R}({\boldsymbol{I}}-{\boldsymbol{S}}){{\boldsymbol{K}}}_{0})}^{-1}({\boldsymbol{S}}+\frac{\alpha }{R}({\boldsymbol{I}}-{\boldsymbol{S}}){{\boldsymbol{K}}}_{0}).$$

(17)

Unlike the single inhomogeneity problem without interfacial thermal resistance, the modified localization tensor B^M should be obtained by decomposing the original problem into three independent problems²². We obtain the relationship between the heat flux ${{\boldsymbol{q}}}^{{\rm{inh}}}$ within a spherical inhomogeneity and the far field heat flux ${{\boldsymbol{q}}}^{{\rm{ext}}}$ as ${{\boldsymbol{q}}}^{{\rm{inh}}}={B}^{M}{{\boldsymbol{q}}}^{{\rm{ext}}}$, where ${{\boldsymbol{B}}}^{M}={[{\boldsymbol{C}}{({{\boldsymbol{C}}}^{{\boldsymbol{M}}})}^{-1}+{\boldsymbol{C}}{{\boldsymbol{K}}}_{0}({{\boldsymbol{K}}}_{1}^{-1}-{{\boldsymbol{K}}}_{0}^{-1})]}^{-1}$. C^M is the modified conjugate Eshelby tensor which is given as ${{\boldsymbol{C}}}^{{\boldsymbol{M}}}=-\,{{\boldsymbol{K}}}_{0}{{\boldsymbol{S}}}^{{\boldsymbol{M}}}{{\boldsymbol{K}}}_{0}^{-1}+{\boldsymbol{I}}$. After substituting the conjugate Eshelby tensor and the modified conjugate Eshelby tensor into the modified localization tensor equation, the heat flux within the inhomogeneity in the presence of interfacial thermal resistance can be obtained as follows,

$${q}_{I}^{{\rm{i}}{\rm{n}}{\rm{h}}}=\frac{\kappa }{{S}_{II}\kappa +{k}_{I}(1-{S}_{II})(1+\alpha \kappa /R)}{q}_{I}^{{\rm{e}}{\rm{x}}{\rm{t}}}$$

(18)

The theoretical predictions of heat fluxes along three directions match very well with the FEA calculation results with the same boundary conditions for a wide range of interfacial thermal resistance $\alpha$, as shown in Fig. 4. At the limit of zero resistance, i.e., $\alpha \to 0$, the heat flux within the inhomogeneity converges to the zero interfacial resistance case depicted in Fig. 3. At the opposite limit, as $\alpha \to \infty$, the heat flux within the inhomogeneity reduces to zero, which implies that the infinite interfacial thermal resistance is equivalent to a void or a perfect thermal insulator as an inhomogeneity. Hence, the effective conductivity approaches that of a porous medium.

Figure 4

Normalized heat flux within the inhomogeneity with respect to the interfacial thermal resistance in a single inhomogeneity problem. The isotropic thermal conductivity of the inhomogeneity is 10 (W/mK).

We then derive a closed form solution for the effective thermal conductivity based on a mean field micromechanics model, the Mori-Tanaka method. Following the previous study²², the effective thermal conductivity of a composite with interfacial thermal resistance can be determined by

$${{\boldsymbol{K}}}_{{\rm{eff}}}={\{{{\boldsymbol{K}}}_{1}^{-1}+\frac{\alpha }{R}{\boldsymbol{I}}-{c}_{0}[{{\boldsymbol{K}}}_{1}^{-1}-{{\boldsymbol{K}}}_{0}^{-1}+\frac{\alpha }{R}{\boldsymbol{I}}]{\boldsymbol{L}}\}}^{-1}$$

(19)

where ${\boldsymbol{L}}={\{{c}_{1}{[{\boldsymbol{C}}{({{\boldsymbol{C}}}^{{\boldsymbol{M}}})}^{-1}+{\boldsymbol{C}}{{\boldsymbol{K}}}_{0}({{\boldsymbol{K}}}_{1}^{-1}-{{\boldsymbol{K}}}_{0}^{-1})]}^{-1}+{c}_{0}{\boldsymbol{I}}\}}^{-1}$ and R is the radius of the inhomogeneity. The effective thermal conductivity of the composite is finally obtained in the closed form as

$${K}_{I}^{{\rm{e}}{\rm{f}}{\rm{f}}}=\frac{{k}_{I}[\kappa +(1-{S}_{II}){c}_{0}((1+\alpha \kappa /R){k}_{I}-\kappa )]}{[{k}_{I}(1+\alpha \kappa /R)-{S}_{II}{c}_{0}((1+\alpha \kappa /R){k}_{I}-\kappa )]}$$

(20)

Although we consider a thermally isotropic inhomogeneity in the final solution for the sake of simplicity, the effective thermal conductivity with anisotropic inhomogeneity can be easily predicted by using an anisotropic K₁.

We plot the effective thermal conductivities, ${K}_{1}^{{\rm{eff}}}$, ${K}_{2}^{{\rm{eff}}}$ and ${K}_{3}^{{\rm{eff}}}$, as a function of the inhomogeneity’s volume fraction ${c}_{1}$ and the interfacial thermal resistance α, as shown in Fig. 5. At the limit of zero interfacial resistance, i.e., ${\rm{\alpha }}\to 0$, Eq. (20) becomes identical to Eq. (12), where we assumed no interfacial resistance. At the opposite limit of ${\rm{\alpha }}\to \infty$, because no heat flux is permitted within the inhomogeneity, the effective conductivity tensor converges to that of a porous medium,

$$\mathop{{\rm{l}}{\rm{i}}{\rm{m}}}\limits_{\alpha \to {\rm{\infty }}}\frac{{k}_{I}[\kappa +(1-{S}_{II}){c}_{0}((1+\alpha \kappa /R){k}_{I}-\kappa )]}{[{k}_{I}(1+\alpha \kappa /R)-{S}_{II}{c}_{0}((1+\alpha \kappa /R){k}_{I}-\kappa )]}=\frac{(1-{S}_{II}){c}_{0}{k}_{I}}{1-{S}_{II}{c}_{0}}$$

(21)

Here, the right-hand side is identical with the zero interfacial resistance solution, Eq. (12), with ${\rm{\kappa }}=0$ (i.e., porous medium solution)

$${\frac{{k}_{I}[\kappa +(1-{S}_{II}){c}_{0}({k}_{I}-\kappa )]}{[{k}_{I}-{S}_{II}{c}_{0}({k}_{I}-\kappa )]}|}_{\kappa =0}=\frac{(1-{S}_{II}){c}_{0}{k}_{I}}{1-{S}_{II}{c}_{0}}$$

(22)

These two limiting cases of ${\rm{\alpha }}\to 0$ and ${\rm{\alpha }}\to \infty$ define the upper bound and the lower bound of the thermal conductivity values, respectively. In plotting Fig. 5, we set ${k}_{1}=1,{k}_{2}=2,{k}_{3}=3\,({\rm{W}}/{\rm{mK}})$, $\kappa =10\,({\rm{W}}/{\rm{mK}})$, and $R=1\,mm$. Because the inhomogeneity is more conductive than the matrix ($\kappa > {k}_{1},{k}_{2},{k}_{3}$), the effective thermal conductivity increases with the volume fraction in the range where the interfacial thermal resistance is low. However, at high enough interfacial thermal resistance, the effective thermal conductivity decreases with the volume fraction because the heat flux through the inhomogeneity is significantly limited. We calculate the critical interfacial thermal resistance that makes the effective thermal conductivity of the composite ${K}_{I}^{{\rm{e}}{\rm{f}}{\rm{f}}}$ identical to the thermal conductivity of the matrix ${k}_{I}$, as follows:

$${\alpha }^{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}}}=R(\frac{1}{\kappa }-\frac{1}{{k}_{I}}).$$

(23)

When the interfacial thermal resistance is higher than the critical value, the effective thermal conductivity of the composite decreases as the volume fraction increases, even though the inhomogeneity is thermally more conductive than the matrix. We note that the critical interfacial resistance decreases for smaller particles, because the interface area increases with the decreasing size of the inhomogeneity for a given volume fraction.

Figure 5

Effective thermal conductivity of composite ((a) ${K}_{11}^{{\rm{eff}}}$, (b) ${K}_{22}^{{\rm{eff}}}$, (c) ${K}_{33}^{{\rm{eff}}}$) having interfacial thermal resistance as a function of volume fraction of inhomogeneity. The thermal conductivities of the matrix and inhomogeneity are ${k}_{1}=1,{k}_{2}=2,{k}_{3}=3$ and $\kappa =10$ (W/mK), and the radius of the particle is 1 (mm).

We validate our analytical solution presented in Eq. (20) by comparing it with the effective thermal conductivity calculated numerically by FEA, as depicted in Fig. 6. We obtain the effective conductivity by averaging the results from 10 independent RVEs, each containing multiple spherical inhomogeneities that are randomly distributed within a cube (see Fig. 7). We assign a uniform temperature boundary condition at two parallel surfaces while applying a periodic boundary conditions along the other two directions. We then compute the heat flux to obtain the effective thermal conductivity of each RVE. As shown in Fig. 6(a–c), the FEA results match very well with our solution for up to 10% of the inhomogeneity volume fraction for a wide range of interfacial thermal resistances. We also investigate the effect of the inhomogeneity’s size at a fixed volume fraction of 5% and an interfacial thermal resistance (10⁻³ m² K/W), as depicted in Fig. 6(d–f). When the interfacial resistance is absent, both theoretical predictions and numerical results find that the effective conductivity is independent of the size of inhogeneity. In contrast, in the presence of interfacial resistance, the effective thermal conductivity decreases as the radius of inhomogeneity decreases, because the interface fraction is bigger for smaller inhomogeneities for a fixed volume.

Figure 6

Effective thermal conductivity ((a)${K}_{1}^{{\rm{eff}}}$, (b)${K}_{2}^{{\rm{eff}}}$, (c)${K}_{3}^{{\rm{eff}}}$) of particle reinforced composite for different interfacial thermal resistances. The radius of the particle used in (a), (b) and (c) is 1 (mm). The effective thermal conductivity ((d)${K}_{1}^{{\rm{eff}}}$, (e)${K}_{2}^{{\rm{eff}}}$, (f)${K}_{3}^{{\rm{eff}}}$) as a function of radius of the particle under fixed volume fraction of 5%. The thermal conductivities of the orthotropic matrix and inclusion are ${k}_{1}=1,{k}_{2}=2,{k}_{3}=3$ and $\kappa =10$ (W/mK) respectively.

Figure 7

Mesh configuration of inhomogeneities in representative volume element for FEA. The volume fraction is 5% and the particle radius is 1 (mm).

Conclusion

In conclusion, we have investigated the heat conduction problem of composites with orthotropic matrices and a spherical inhomogeneities in the presence of interfacial thermal resistance. We derive the modified Eshelby tensor of the eigen-intensity problem as well as the effective thermal conductivity based on a micromechanics approach by considering the interfacial thermal resistance effect, and validate our solution against FEA calculation results. We also demonstrate that the effective conductivity solution has the correct limiting behaviour at both the zero and infinite interfacial thermal resistance limit. The solution in the present paper is applicable to the composites with either transversely isotropic or isotropic matrices and an inhomogeneity with an arbitrary thermal conductivity tensor. We plan to extend the present study by considering the size effects of nanoscale inclusions for nanocomposites⁴⁹ in obtaining an analytic solution and by coupling molecular dynamics simulations of ceramic composites or polymer composites with the analytic solution. We believe that our study can provide an effective and accurate way of predicting the thermal conduction of composites, and it can be applied to better design technologically important materials such as polymer-based composite and thermoelectric materials.